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ABSTRACT  ' 

This  report  provides  a  brief  Introduction  to  the  use  of  MACSYMA,  a 
symbolic  manipulation  program  for  mathematics. 

The  first  Chapter  outlines  invoking  Maesyma,  running  tutorials  and  demos, 
saving  transcripts  and  programs,  editing  and  plotting  facilities. 

The  second  Chapter  provides  12  examples  of  Maesyma  use.  The  following 
topics  are  touched  upon:  factorization,  integration,  matrix  multiplications, 
eigenvalue  problem,  infinite  series,  recursion  definitions  and  relations, 
solving  systems  of  polynomial  equations  and  ODEs. 

Some  excunples  in  Chapter  2  and  the  intricate  exan^le  given  in  Chapter  3 
demonstrate  the  possibilities  of  programming  in  Maesyma.  - 


AMS  (MOS)  Subject  Classification:  68Q40 

Key  Words:  maesyma,  symbolic  manipulation,  symbolic  programming 


Supported  in  part  by  the  U.  S.  Army  Research  Office  under  Grant  No.  DAAL03-87- 
K-0028  and  the  Air  Force  Office  of  Scientific  Research  under  Grant  No.  85- 
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MACSYMA  AT  CMS.  VERSION  309.3 


W.  Hereman,  Y.  Nagel  and  J.  strlkwerda 


Chapter  1 

Introduction 


Macsyma  is  a  symbolic  manipulation  language  for  mathematics  written  in  Lisp.  It  was  originally 
developed  by  the  Mathlab  group  of  the  MIT  Computer  Science  laboratory.  It  is  now  supported  and 
distributed  by  Symbolics,  Inc.  of  Cambridge,  Massachusettes.  It  can  perform  complicated  math¬ 
ematical  computations  involving  integration,  differentiation  and  matrix  or  function  manipulations. 
It  can  also  plot  curves  on  suitable  graphics  terminals. 

This  document  is  an  introductory  tutorial  to  using  Macsyma  on  the  CMS  Vax  11/780.  It  is  not 
intended  to  be  a  comprehensive  guide  or  reference.  Basic  references  for  Macsyma  are: 


Vax  Unix  Macsyma 
Reference  Manual 
Version  11 
Symbolics,  Inc. 
Cambridge,  Mass 


Applications  of 
Macsyma  to  Calculations 
in  Applied  Mathematics 
M. Hussain  k  B.  Noble 
GE  Report#83CRD054 


Macsyma:  A  Program  for 
Computer  Algebraic 
Manipulation 

Naval  Underwater  Systems 
Center,  Newport,  R.l. 


The  first  book  is  the  basic  guide  to  Macsyma.  It  is  arailable  at  the  MACC  documentation  center. 
It  is  not  a  good  tutorial  although  it  does  contain  a  list  of  other  references.  The  other  two  books 
have  many  examples  suitable  for  learning  Macsyma.  Some  of  these  examples  are  quite  complicated. 

1.1  Invoking  Macsyma 

To  use  Macsvma  on  the  CMS  VAX  you  must  first  log  in  to  the  computer  (see  the  System  Manager 
if  you  do  not  have  an  account).  When  you  first  login  you  will  see  the  VMS  t  prompt.  The  CMS  VAX 
has  Eunice,  a  unix  emulator,  as  well  as  the  basic  V'MS  operating  system.  You  must  run  Macsvma  from 
Eunice.  To  enter  Eunice,  you  first  type  unix  or  eunice.  You  will  then  get  the  unix  prompt  -  usually 
a  %  or  — >.  Now  type  in  /usr/maesyma  to  start  macsyma.  Macsyma  will  return  the  (cl)  prompt. 

I  S  Unix  ^1 

I  %  /usr/maesyma  | 

i  (c»)-  ~  i 

As  you  work  through  an  example,  Macsyma  prints  (cn).  where  “n”  is  an  integer,  for  user  input 
lines  and  (dn)  or  (en)  for  its  own  replies.  The  (en)  lines  denote  intermediate  calculations  used 
by  Macsyma  for  especially  long  computations  or  lists  of  output.  Lines  are  numbered  consecutively 
and  can  be  used  to  recall  previous  commands  or  responses.  All  Macsyma  commands  must  end  with 
either  a  I  or  a  ;.  Use  the  f  if  you  don’t  want  a  d-Iine  to  echo  your  input  or  to  avoid  intermediate 
output.  Lines  enclosed  in  double  quotes  or  /*  ■  •  •  •/  are  treated  as  comments  (see  Example  12 
in  Chapter  2).  You  can  include  such  lines  (terminate  them  with  a  I)  in  your  Macsyma  session  to 

Supported  in  part  by  the  U.  S.  Army  Research  Office  under  Grant  No. 
DAAL03-87-K-0028  and  the  Air  Force  Office  of  Scientific  Research  under 
Grant  No.  85-0263. 
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remind  you  of  the  purpose  of  the  calculation.  See  the  next  section  for  instructions  on  saving  your 
output  in  a  transcript  file.  The  percent  sign,  %.  can  be  used  to  refer  to  the  last  expression.  In  the 
examples,  boldface  indicates  the  characters  typed  in  by  the  user.  Normal  letters  show  characters 
produced  by  the  computer. 

Some  special  characters  are; 


for  the  imaginary  number, 

%e 

for  the  base  of  the  natural  logarithms, 

%pi 

for  X, 

% 

refers  to  expression  on  previous  line 

%th(n) 

refers  to  expression  on  n‘^  previous  line 

inf 

for  oo. 

The  Macsyma  part  command  allows  you  to  refer  to  a  part  of  the  n‘^  previous  expression  (see 
the  examples  in  Chapter  2.) 

Note  that  the  exponential  function  with  argument  x  can  be  typed  in  two  diiferent  ways 

%  ex 
or  exp(x) 

To  exit  from  Macsyma  use  quit();.  Ctrl  c  will  force  an  interrupt.  You  can  then  type  h  to  see 
what  options  you  have.  The  most  useful  ones  are  reset  which  returns  Macsyma  to  a  c-line  and 
quit  which  will  make  you  exit  from  Macsyma  and  return  to  Eunice.  Use  ctrl  c  or  Ctrl  y  to  end  a 
runaway  calculation.  Ctrl  c  will  toggle  you  into  or  out  of  LISP. 


1.2  Tutorial  and  Demos 

An  on-line  tutorial  for  Macsyma  is  available.  To  access  it,  type  in  /usr/macsyma  then 

(cl)priineT();[ 

A  menu  will  appear.  Typing  in  the  number  of  an  item  on  the  menu  will  cause  a  Macsyma  script  to 
appear  with  appropriate  examples  and  instructions. 

Macsyma  also  has  two  built-in  libraries  of  examples.  The  share  library  /usr/mac309/share 
can  be  used  with  the  example  command.  If  at  any  point  (after  typing  in  /usr/macsyma  you 
would  like  to  see  how  a  function  can  be  used,  type  in  example(/unclion)  for  a  brief  tutorial.  For 
instance, 

(cn)  example(ode2); 

will  yield  a  brief  introduction  (with  examples)  of  the  odc2  solver.  To  read  a  description  of  a  command, 
option,  or  concept  type,  e.g. 

(cn)  deBcribe(stringout); 

will  give  you  a  description  of  the  command  etrtngout.  See  the  list  of  special  forms  in  the  manual 
index  for  other  possible  commands.  This  is  a  list  of  the  .dem  files  in  /u$r/mac309/share; 


absimp 

dimen 

nchrpi 

rncomb 

antid 

eigen 

ode2 

sqdnst 

dblint 

facexp 

optmii 

submac 

delta 

fourie 

pfaff 

trgsmp 

differ 

Irats 

polsol 

vect 

The  demonstration  library  /usr/macSOQ/demo  also  contains  many  examples.  You  can  select 
demos  from  this  library  by  typing 

(cl  )batch(”/usr /mac309/deino/ /uncttDn.dem”; 

or 

(cl)demo(”/usr /mac3O0 /demo//u7ictioTi.dem”  )7j 

The  batch  command  causes  Macsyma  to  run  through  the  demo  without  stopping  (as  though 
you  were  looking  at  a  movie).  With  the  demo  command.  Macsyma  will  pause  after  each  display 
and  show  an  underscore  _  cursor.  Hitting  the  Return  key  will  cause  it  to  continue.  To  get  a  list 
of  the  files  in  these  two  directories  type  Is  /usr/mac309/demo  or  Is  /usr/mac309/share  from 
Unix  -  i.e.  at  the  %  prompt. 

Here  is  a  list  of  the  .dem  files  in  /usr/mac309/demo: 


algfac 

differ 

matrix 

series 

aJgsys 

difficult 

mpg 

short 

array 

ezgcd 

nalgfc 

simpl 

ball 

factor 

newfac 

sin 

begin 

gen 

nisimp 

solder 

e2cyl 

iap 

pfacto 

solve 

cf 

int 

plot 

specfh 

combin 

laplac 

poplot 

subscr 

cyl2c 

legen 

qpr 

sum 

cyl2s 

limit 

radcan 

taylor 

decomp 

macex 

ratsimp 

trace 

defint 

macro 

ratsub 

trig 

demall 

maneg2 

risch 

demtest 

maneg3 

rough 

dice 

mat 

rpart 

1.3  Saving  Macsyma  Output 

You  can  save  your  Macsyma  output  as  either  a  TRANSCRIPT  file  or  a  LISP  file.  The  TRANSCRIPT 
file  contains  all  the  c-  and  d-lines  which  appear  on  your  terminal  during  an  interactive  session  with 
Macsyma.  The  LISP  file  contains  LISP  commands  generated  by  Macsyma  and  can  be  used  to 
playback  or  restart  your  last  interactive  session 

1.3.1  Transcript  file 

To  save  just  the  screen  output  (w'ith  displayed  equations)  do  the  following,  type  writefile  at  the 
beginning  and  closefile  at  the  end  of  a  computation  as  in  the  following  example 

(cm)  writefile(foo);  -  logs  to  file  ’’foo” 


(cn)  closefile();  -  closes  log  file  "foo” 

Get  out  of  macsyma.  Send  "foo”  to  the  printer,  "foo”  cannot  be  reloaded  to  Macsyma.  See 
section  10-17  of  the  Macsyma  manual.  If  you  type  quit();  without  doing  a  closefile();,  a  transcript 
file  will  be  created  and  saved  in  your  directory  anyway  If  you  forget  to  create  a  writefile  when  you 
begin  your  macsyma  session,  you  can  create  one  at  any  time. 

3 


(cm)  writefile(foo); 
(cm+1)  playback(all): 

(cn)  closefileO; 


-  logs  to  file  "foo” 

-  repeats  your  entire  session 

-  closes  log  file  ”foo” 


If  we  type  for  line  (cm-^1) 

'  playback(c35,c39); 

Macsyma  would  only  playback  lines  c35  thru  c39.  This  allows  you  to  be  more  selective  about  what 
you  save. 

1.3.2  LISP  file 

To  save  LISP  commands,  suitable  for  reloading  or  playing  back  your  session  see  section  10-18  of  the 
Macsyma  manual. 

Enter  various  Macsyma  commands.  At  end  of  session  do  this: 


(cm)  Bave(soo,all);  -  saves  all  Macsyma  commands  in  Lisp 

in  the  file  ”soo” 

(cn)  quit();  Exit  from  Macsyma 

If  you  qaitO;  before  typing  8ave(soo);,  you  will  not  be  able  to  recover  the  LISP  transcript  of 
your  session.  To  reload  the  LISP  file  just  re-enter  Macsyma  and  issue  this  as  the  first  command: 


(cl)  loadfile(8oo);  -  Reloads  "soo”  Lisp  file  into  Macsyma 

(c2)  play  back  (soo);  -  Recreates  your  last  session  with  ”soo” 

This  is  optional.  It’s  not  required  to  do 
further  computations.  See  10-4  in  manual 


An  efficient  way  to  save  a  single  line  of  output  (say  the  final  result  of  your  macsyma  computations) 
in  LISP  form  is  shown  in  the  following  example: 


(cl)  depeiids([f,v].[xl,x2,x3]): 

(c2)  f:(xl)'2*(x2)*(x3)  +  (xl)*(x2)-2*(x3)  - 
(c3)  v(f):[difF(f,xl),difr(f,x2),diff(f.x3)); 

(d3)  •  •  • 

(c4)  stringout(grad.c4): 

(c5)  quit(); 

%  /usr/macsyma 

(cl)  f:xl  -f  x2'2  -i-  x3*3; 

(c2)  batch(grad): 

(c3)  v(f):[difr(f.xl),difr(f,x2),difr(f,x3)): 


(xl)*(x2)»(x3)-2; 

-  v(f)  is  the  gradient  of  f 

-  Macsyma’s  answer  is  on  this  line 
save  (c4)  in  Macsyma  file  grad 
out  of  macsyma 

reenter  macsyma 
declare  f 

input  definition  of  grad 
machine  returns  answer 


You  can  use  the  sln’n^out  command  to  save  selected  c-,  d-,  or  e=  lines.  The  most  useful  possibilities 

are 


stringout(/i/enamc.ck) 

5tringout{filcnam(.cj,ck) 

stringout(^lcnamf,[j.k|) 


saves  line  (ck) 

saves  lines  (cj)  1'  (ck) 

saves  lines  (cj)  thru  (ck) 


You  can  examine  the  contents  of  grad  by  tvping  type  grad  at  the  VMS  dollar  sign  prompt.  You 
can  also  import  grad  into  Macsyma  via  the  command  batch(grad); 


k***.t*- 


'Id 


1.3.3  Combining  two  transcripts 

Suppose  we  have  saved  only  the  Lisp  commands  in  the  file  ”soo”.  Now  we  want  to  get  the  displayed 
equations.  Do  this 

(cl)  writeftle(doo): 

(c2)  loadfile(soo); 

(c3)  pla)'back():  don’t  need  to  mention  soo 

(c4)  closefileO; 

(c5)  quit();  get  out  of  Macsyma 

%  Ipr  doo  send  “doo”  to  laser  printer 

1.4  Running  Jobs  in  the  Background 

MACSYMA  version  309.3.  unlike  the  older  version  of  MACSYMA  cannot  be  called  from  VMS.  This 
means  that  macsyma  cannot  be  run  on  the  VMS  batch  queues.  Instead,  you  can  run  a  macsyma 
job  in  the  background  as  follows: 

1.  First,  use  an  editor  to  create  a  file,  called,  say,  JUNK,  which  contains  the  MACSYMA  com¬ 
mands  you  would  have  typed  at  the  terminal  if  you  were  using  MACSYMA  interactively.  A 
sample  file  might  consist  of  these  lines 


(cl)  writeiile(trash)t 
(c2)  f(x);=  l/tan(x); 
(c3)  difF(f(x),x); 

(c4)  closefUeO; 

(c5)  quitP; _ 


2.  Notice  that  JUNK  and  TRASH  are  different  files.  JUNK  is  a  file  you  create  with  an  editor. 
TRASH  is  the  file  Macsyma  will  create  to  save  the  output  of  your  program.  There  are  three 
possible  ways  to  pipe  JUNK  into  MACSYMA 

(a)  /usr /macsyma  <  junk 

This  is  closest  to  running  the  job  interactively. 

(b)  /usr/macsyma  <  junk.  & 

This  pipes  JUNK  into  MACSYMA  and  puts  it  into  the  background 

(c)  nice  4-  <PN>  /usr/macsyma  <  junk.  & 

This  pipes  JUNK  into  MACSYMA,  puts  it  into  the  background, 
and  assigns  it  PRIORITY  PN.  Acceptable  choices  for  PN  are 
integers  in  the  range  10  thru  20.  The  higher  the  number, 
the  lower  the  priority;  10  is  the  default. 


1.5  Using  the  “vi”  Editor  with  Macsyma 

The  default  editor  for  MACSYMA  is  TECO.  If  you  wish  to  use  a  different  editor,  say  “vi”  with 
Macsyma,  create  a  file  called  MACSYMA. MAC  (with  edt  or  vi  or  any  editor  you  like)  which  contains 
the  following  lines: 

/•-•-Macsyma-*-*/ 
ioad(‘‘ucb//newedit.o”); 
editorcom:“vI  macsyma. buT'; 


!««  «•« 


lJ  'Ll'kl'  •.•■  •.t'a^  ‘a 


MACSYMA.MAC  is  an  initialization  file  for  MACSYMA.  It  must  be  in  your  home  directory  (the 
directory  you  are  in  when  you  first  log  on).  Next,  type  the  following  V'MS  command 


I  S  set  term/noline-edit 


Follow  the  procedure  given  above  for  invoking  MACSYMA  interactively.  Now  suppose,  that 
while  I  am  in  macsyma,  I  type  the  following  at  command  line  (cl): 


Ij;  -  no  <RET> 


If  I  now  hit  the  <ESC>  and  <RET>  keys,  1  will  put  this  line  into  the  editor,  suitable  for  editing 
with  vi.  If  by  mistake,  1  enter  LISP,  I  can  type: 


to  leave  MACSYMA  altogether  or  Ctrl  s  to  return  from  LISP  to  (cn). 

If  I  have  gone  considerably  beyond  line  (cl)  but  want  to  re-edit  it  to  use  it  again,  I  would  type: 


string(cl);  <ESC>  <RET> 


and  line  (cl)  will  be  put  back  into  the  buffer.  Although  this  technique  will  allow  you  to  edit  long 
lines,  it  is  better  to  enter  long  lines  by  breaking  them  up  into  shorter  pieces  and  combining  them 
step  by  step. 


1.6  Plotting  on  a  Tektronix  with  Macsyma 


Macsyma  will  plot  on  any  terminal  which  emulates  a  Tektronix  4010  or  4014  terminal.  At  CMS 
some  of  the  V'isual  102  graphics  terminals  do  this.  IBM  PCs  with  suitable  graphics  cards,  graphics 
monitors  and  Tektronix  emulation  software  can  do  this  also.  You  must  first  put  your  terminal  or 
PC  in  Tektronix  emulation  mode  before  attempting  to  plot.  Here  is  an  example  using  the  Macsyma 
function  plot2 


First,  put  your  PC  or  terminal  into  Tektronix  emulation  mode 
S  Unix 

%  set  term=4014 
%  /usr/macsyma 
(cl)  plotmode:tek; 

(c2)  plot2([x-fl,x'2-t-l].x,-l,l); 


i; 


Chapter  2 

Examples 


You  can  try  the  examples  (below)  to  get  an  easy  introduction  to  Macsyma 

2.1  Example  1  -  polynomial  factoring 

To  factor  the  polynomial  x®®  -t-  1  enter: 

(cl)  x*99  +  1; 

(c2)  factor(%); 


2.2  Example  2  -  repeating  a  command 

(cl)  f(x);=  l/tan(x); 

(c2)  uitegrate(f(x),x); 

(c3)  f(x):=  x/((x+l)*(x'2+l)*2); 

(c4)  ’  ’(c2); 

(c5)  quit(); 

Note  the  use  of  the  two  single  quotes  on  line  (c4)  to  repeat  the  command  on  line  (c2). 


2.3  Example  3  -  matrix  multiplication 

Macsyma  has  two  different  kinds  of  matrix  products.  An  asterisk  *  will  produce  an  element-wise 
product,  while  a  period  •  yields  a  true  matrix  product.  Here  is  an  example 
(cl)  enteriTiatrix(2,2); 

(c2)  M:%; 

Macsyma  will  assign  the  label  m  (or  M  -  Macsyma  is  case-insensitive  to  the  matrix  above 
(c3)  M*M; 

Macsyma  will  square  each  element  of  the  matrix 
(c4)  M**2; 

This  will  yield  the  same  result  as  line  (c3) 

(c5)  M.M; 

Macsyma  will  return  the  true  matrix  product 
(c6)  M"2; 

This  will  yield  the  same  result  as  line  (c5) 


2.4  Example  4  -  matrix  eigenvalues 


To  find  the  eigenvalues  of  the  matrix 

a  “) 

V 

enter: 

(cl)  entermatrix(2,2); 

Is  the  matrix  1. Diagonal  2. Symmetric  3, Antisymmetric  4. General 
Answer  1,  2,  3,  or  4 

4; 

Row  1  Column  1:  1; 

•  •  •  Put  in  remaining  entries.  Remember  the  semicolon! 

(c2)  eigenvalues(%); 
will  give  you  the  eigen\'alues. 


2.5  Example  5  -  infinite  sums 

To  evaluate  the  sum  o 

(cl)  assuine(  abs(x)  <  1); 

(c2)  sum(x'j  ,  j,  0,  inf); 

(c3)  simpsumttruet 
(c4)  ev(d2,suin); 

(c5)  forget(  abs(x)  <1)1 
(c6)  ev(d2,sum): 

Is  abs(x)  -  1  positive,  negative,  or  lero? 
negative; 


1 

1  -  X 

(c7)  ev(d2,sum); 

Is  Bbs(x)  -  1  positive,  negative,  or  zero? 
aero; 

undefined 

(c8)  ev(d2,sum); 

Is  abs(x)  -  1  positive,  negative,  or  zero? 
positive; 

inf 

If  you  don’t  put  in  the  assume  instruction  (line  (cl))  you  will  be  queried  on  the  size  of  z  as  in 
lines  (c6)  thru  (c8). 


2.6  Example  6  -  adding  and  subtracting  series 

This  example  shows  how  to  add  and  subtract  infinite  sums. 

(cl)  sum(a[k]  +  b[k].k.0,inf); 

(c2)  suin(b[k].k.0anf); 

(c3)  dl  -  d2; 

(c4)  timpsunitfalse; 

(c5)  intosum(%); 


(c6)  suincontract(%); 

The  aimpsum  option  must  be  set  to  false  in  order  for  intosum  to  work  properly.  Without  the 
intosum  command  the  sumconiTaci  will  not  do  anything,  because  of  the  minus  sign  in  front  of  the 
second  sum.  If  the  sums  were  added  then  there  is  no  need  for  the  intosum.  For  example, 

(cl)  sum(a[k]  -  b[k],k,0,inf); 

(c2)  suin(b[k],k,0,inf); 

(c3)  dl  +  d2; 

(c4)  sumcontract(%); 

will  produce  the  same  answer  as  above,  i.e.  the  sum  of  the  o*. 


2.7  Example  7  -  a  bug  in  macsyma 


To  evaluate  the  integral 


Jo  +  a 


(cl)  ’iiitegrate(l/(x'2  4-  a),  x,0,inf); 

Macsyma  will  respond  by  displaying  the  integral. 

(c2)  ev(%^ntegrate,a=l); 
will  produce  the  result: 

(d2)  %pi  /2 

The  use  of  ’integrate  produces  the  expression  of  the  integral,  while  using  integrate  would  cause 
evaluation  of  the  expression.  Continuing  this  example,  we  may  wish  to  evaluate  the  derivative  of 
(cl)  with  respect  to  a. 

(c3)  difr(cl^); 

(c4)  assuine(a>0); 

(c5)  ev(d3:integrnte); 

Macsyma  gives  a  strange  result  for  this, 


If  we  now  type 

(c6)  lognegint:true9 
Macsvma  will  return 


which  hats  the  wrong  sign 


2.8  Example  8  -  checking  an  equation 

This  example  shows  lhal  the  polynomials,  Pfc(j).  dcf.ned  by 

00 

exp(y^)  =  y^/p*(x) 

*=o 

satisfy  the  equation 

xPk(x)  -  kpk{x)  -  {k  -  l)p*-,(z). 
We  will  only  verify  it  for  k  from  0  through  6 


(cl)  exp((x*s)/(l-s)); 

(c2)  taylor(%,s,0,6); 

(c3)  expand(%); 

(c4)  for  i  from  0  thru  6  do  a[i]:ratcoefr(d3,s4); 
note  that  if  we  had  written  instead 

for  i  thru  6  do  a[i]:ratcoefF(c3.s,i); 
macsyma  would  accept  this  command  but  do  only  terms  1  thru  6 
(c5)  for  i  from  0  thru  6  do  b[l]:x*difr(a[i],x)  -  i*a|i]  +(i-l)*a[i-l]; 

At  this  point  the  command 
(c6)  for  i  from  0  thru  6  do  ldisplay(b[i]); 
will  display  polynomials  but  they  will  not  be  simplified  yet.  Use 
(c7)  for  i  from  0  thru  6  do  b[i]:expand(b[i]); 

(c8)  for  i  from  0  thru  6  do  Idisplay(b[i]); 

to  see  all  the  expected  zeroes.  Note  that  the  command 

(c7)  for  i  from  0  thru  6  do  expand(b[i]); 

would  do  no  good  since  the  result  of  the  expand  is  not  kept  anywhere.  Instead  of  expand  one  could 
also  use  ratsimp  or  gfactor 

2.9  Example  9  -  solving  systems  of  equations 

The  solve  command  can  be  used  on  systems  of  linear  or  polynomial  equations  in  several  variables. 
To  find  all  complex  solutions  to  x  +  y  +  z  =  2,  x*  +  y*  +  z*  =  26,  +  y®  -t-  z®  =  38,  enter: 

(cl)  x+y+*=2; 

(c2)  x**2+y**2+***2=26; 

(c3)  x**3+y**3+***3=38; 

(c4)  •olve((dl,d2,d3],[x,y.i]); 

(d4)  [[z  =  -  3,  y  =  1,  z  =  4] ,  [x  =  -  3,  y  =  4,  z  =  l] , 

[x  =  4,y  =  -3,  z  =  l],  [x  =  4,y  =  l,z=-3],  [x  =  l,y  =  -3,z  =  4], 

[x  =  1,  y  =  4,  z  =  -  3]] 

The  brackets  indicate  the  solutions  ate  stored  in  an  array.  To  remove  the  first  solution  enter: 

(c5)  d4[l]i 

(dS)  [x  =  -  3,  y  =  1,  I  =  4] 

(c6)  d5[l]; 

(d6)  X  =  -  3 

(c7)  d5l2]; 

(d7)  y  =  1 

(c8)  d5[3); 

(d8)  f  =  4 

To  obtain  parts  of  an  equation  or  expression  use  the  ‘part'  command. 

For  example,  enter 
(c9)  part(d6.1); 

(d9)  X 

(clO)  part(d6.2); 

(dlO)  -  3 

This  example  is  due  to  Paul  Terwilliger 


10 


2.10  Example  10  -  Pascal’s  triangle  and  recursive  defini¬ 
tions 

In  this  example  we  generate  Pascal's  Triangle. 

Define  the  binomial  coefficient  recursively  by 

(cl)  c(n,k):=  if  (k  =  -1  or  k=  n  +  1)  then  0  else  if 
(k  =  0  and  n  =  0)  then  1  else  c(n-l.k-l)+c(n-l,k); 

(dl)  c(n,  k)  :=  IF  k  =  -  1  OR  k  =  n  -  1  THEN'  0 

ELSE  (IF  k  =  0  AND  n  =  0  THEN  1  ELSE  c(n  -  1,  k  -  1)  +  c(n  -  1,  k)) 

We  may  now  evaluate  any  entry  in  the  Triangle,  for  example 

(c2)  c(6,3);  /*  c(6,3)  refers  to  the  fourth  element  of  the  seventh  row  */ 

(d2)  20 

To  print  out  the  first  7  lines  of  the  Triangle,  enter 

(c3)  for  i  :  0  thru  6  do  (array(line,i),for  J:0  thru  i  do  line[i]:c(ij), 
disp(listarray(line))); 

/usr/mac309/maxsrc/arrav.o  being  loaded. 

[1] 

(1,  1] 

U,  2.  1] 

[1,  3,  3,  1] 

[1,  4,  6,  4,  1] 

[1,  5,  10,  10,  5,  1] 

[1,  6,  15,  20,  15,  6,  1] 

(d3)  done 

As  the  variable  i  above  ranged  from  0  to  6,  Maesyma  ran  through  a  procedure  that  printed  out 
the  ith  line  of  the  Triangle.  First,  an  array  ‘line’  was  defined,  with  entries  indexed  by  0,l,...,i.  Then 
as  the  variable  j  ranges  over  these  indices,  linejj]  gets  assigned  the  value  c(i,  j).  After  the  line  is  filled 
its  entries  are  displayed.  This  example  is  due  to  Paul  Terwilliger 


2.11  Example  11  -  solving  ODEs 


This  example  shows  how  Maesyma  handles  ODEs 
Let  us  first  try  a  second  order  ODE 

(cl)  diffeq:  x*2*’diff(y,x,2)-3*x*’difF(y,x)+4*y  =  0; 
(dl) 


(c2)  ode2(difreq,y,x): 

Batching  the  file  /usr/mac309/share/ode2.mac  ... 

Batching  done 
(d2) 

y  =  x^{%k2log(z)  ~  %kl) 


Now  let  us  try  a  third  order  ODE 

(c3)  difTeq:  x*3*’diff(y,x,3)-3*x*’difF(y,x)+S*y  =  0; 
(d3) 


a*—  -  3*—  +  3y  =  0 

di»  dz  ^ 


(c4)  ode2(difFeq,y,x); 


dx^  dx 

Not  a  proper  difTerential  equation 
(d4)  false 

V 

Now  let  us  try  a  first  order  ODE 

(c5)  x-2*'’diff(y,x)  +  3*x*y  =  8m(x)/x; 

(d5) 


(c6)  so]nl:ode2(%,y,x); 
(d6) 


Now  we  give  the  initial  conditions: 
Let  us  require  that  for  x=^ir,  y  =  0 
(c7)  icl(8olnl,x=%pi,y=0); 

(d7) 


Here  is  another  example; 

(c8)  ’diHl(y,x,2)  +  y*'diff(y,x)*3  =  0; 


ax  X 


»  = 


_  %c  —  cos(x) 


y  = 


_  cos(x)  +  I 


(d8) 


— +  y(^)’ 

dx>  ’''dx^ 


=  0 


(c9)  soln2:ode2(%,y,x); 
(d9) 


+  6%fcly 
6 


=  X  +  %k2 


(clO)  rat8iinp(ic2(soln2,x=0,y=0,’difr(y,x)=2)); 
(dlO) 

2y-  3y  _ 


(cll)  bc2(8oln2,x=0,y=l,x=l,y=3); 

(dll) 

y®  -  lOy  3 

— :: - =  *  -  r 


2.12  Example  12  -  programming  in  Macsyma 

This  example  shows  how  to  use  a  file  of  commands  to  execute  the  commands  from  within  macsyma. 
We  will  illustrate  this  by  using  the  algorithm  for  testing  if  a  polynomial  is  a  Schur  polynomial.  A 


Schxir  polynomial  is  one  which  has  all  of  its  roots  inside  the  unit  circle.  An  algorithm  to  test  whether 
or  not  a  polynomial  is  a  Sckur  polynomial  is  as  follows.  Given  a  polynomial  ^n(^)  of  degree  n, 


we  define  ^„(2)  as 


n 

in{i)  = 

t=0 


where  the  bar  on  the  coefficients  denotes  the  complex  conjugate.  The  polynomial  is  defined 

by 

,  ,  ,  «n(0)^„(2)  -  <^n(0)<^„(2) 

=  - • 


The  algorithm  is  based  on  the  fact  that  ^n(^)  i*  »  Schur  polynomial  if  and  only  if  |^„(0)|  is  less  than 
|^„(0)|  and  ^n-i(2)  is  a  Schur  polynomial. 

The  following  macsyma  instructions  have  been  placed  in  a  file  called  sckur.  Comments  are  given 
between  the  delimiters  /*  and  */. 

/*  Code  for  testing  for  Schur  polynomials 
the  polynomial  must  be  p 
the  value  of  n  must  be  the  degree  of  p 
If  aU  the  parameters  t[j]  are  nonnegative  then 
p  is  a  Schur  polynomial  */ 

/•  Run  by  using  batch(schur)  •/ 

schur:  true;  /*  Initialise  schur  to  true  */ 

while  n  >  0  and  schur  =  true  do  ( 
display(n), 

for  i  from  0  thru  n  do  aa[i]:ratcoefl(p,z4)> 

/*  Compute  the  check  on  the  product  of  roots.  */ 

xO:  realpart(aa[0]), 

yO:  imagpart(aa[0]), 

xn:  realpart(aa[n]), 

yn:  imagpart(aa[nl), 

t[n]:  xn"  2  +  yn'  2  -xO'  2  -yO'  2, 
ldisplay(t[n]). 

if  t[n]  <=  0  then  schur:  false  , 

r:  sum{aa[n-i]*z'  idiO,n),  /*  r  is  the  “reverse”  polynomial  of  p  * / 

q:  subst(-%i,%i.r),  /•  q  is  the  tilde  of  p  */ 

bb[0]:  ratcoef(q,z,0), 

r:  bb[0]*p  -  aa[0]*q. 

r:  ratsimp(r), 

p:  ratsimp(r/s). 

Idisplay(p), 
n:  n>l 

);  /•  end  of  while  loop  */ 

Idisplay  ( s  chur  ) ; 

This  set  of  commands  illustrates  the  use  of  the  ‘while  loop’  and  the  ‘for  loops.’  Note  that  the 
commands  within  the  loops  end  with  a  comma. 

As  an  example  of  how  to  use  this  file  we  consider  the  three  polynomials 

(7  +  2i)i*  -  (8  -  0*  +  1. 
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222*  -  202*  -  32  +  1, 
and 

(22  +  *)2*- (20-»)2*-32+ 1. 

We  enter  the  polynomial  and  its  degree  and  then  call  the  file  as  follows. 

(cl)  p:  (r+2*%i)*2-2  -{8.%i)**  +1; 

(c2)  n:  2$ 

(c3)  batch(schur); 

The  use  of  the  $  at  the  end  of  (c2)  supresses  the  display  of  the  value  of  n  on  the  screen.  The 
command  in  (c3)  causes  the  file  sckur  to  be  read  and  the  executed.  The  file  is  also  listed  before  it 
is  executed.  The  final  value  of  the  logical  variable  schur  will  be  displayed.  (In  this  case  it  is  true.) 
Next  we  enter  the  second  polynomial. 

(cll)p:  22*z‘3 -20*2*2  -3*2  +1; 

(cl2)  n:  3S 
(cl3)  batch(schur): 

Again  the  commands  in  the  file  will  be  executed,  and  in  this  case  the  polynomial  is  not  a  Schur 
polynomial.  (  1  is  a  root.)  Finally,  we  enter  the  third  polynomial. 

(c21)  p:  (22  -|-%i)*2*3  -(20-%i)*2*2  -3*2  -fl; 

(c22)  n:  3$ 

(c23)  batch(schur): 

Again  the  commands  in  the  file  will  be  executed,  and  in  this  case  the  polynomial  is  a  Schur 
polynomial. 


Chapter  3 

Intricate  Example 


The  Painleve  Property 

Ref:  J.  Weiss,  J.  Math  Phys.,  24(6),  June  1983,  1405 


3.1  Introduction 

We  say  that  a  partial  differential  equation  has  the  Painleve  property  when  the  solutions  of  the 
pde  are  “single-valued”  about  the  movable,  singularity  manifold.  To  be  precise,  if  the  singularity 
manifold  is  determined  by 

9(ij,...,z«)  =  0  (3.1) 

and  /  =  f(zi,—,Zn)  is  a  solution  of  the  partial  differential  equation,  then  we  assume  that 

/  =  g-f;u,y,  (3.2) 

,=0 

where 

y  =  ff(2i,...,Zn),  tt,  =  tt^(zi,...,r„),  tto  0. 

are  analytic  functions  of  (zj )  in  a  neighborhood  of  the  manifold  (3. 1 )  and  a  is  an  integer.  Substitution 
of  (3.2)  into  the  partial  differential  equation  determines  the  value(s)  of  a  and  defines  the  recursion 

relations  for  Uj.j  =  0,1,2 .  When  the  ansati  (3.2)  is  correct,  the  pde  is  said  to  possess  the 

Painleve  property  and  is  conjectured  to  be  integrable. 


3.2  The  Korteweg-de  Vries  Equation 

The  KdV  equation 

/.  +  //.  ^  =0.  6  e  IR 

possesses  the  Painleve  property.  The  expansion  about  the  singular  manifold  has  the  form 

>=o 


It  is  found  that  the  “resonances”  occur  at 


J  =  -1,4,6. 


(3.3) 


(3.4) 


(3.5) 


IS 


JaV.' 


Resonances  are  those  values  of  j  at  which  it  is  possible  to  introduce  arbitrary  functions  into  the 
expansion  (3.4).  For  each  nontrivial  resonance  there  occurs  a  compatibility  condition  that  must  be 
satisfied  if  the  solution  is  to  have  a  single-valued  expansion.  The  resonance  at  j  =  —  1  corresponds 
to  the  “arbitrary”  function  g  defining  the  singular  manifold  (3.3).  The  compatability  conditions  at 
j  =  4  and  6  are  satisfied  identically. 

From  the  recursion  relations,  we  find: 


j  =  0:  Uo=-12bs|;  (3.6) 

i  =  1  :  It,  =  126s„;  (3.7) 

3  =  2  :  +  46s* 9***  -  36s’,  =  0;  (3.8) 

i  =  3  :  s*«  +  -  s*“3  -  6s****  =  0;  (3.9) 

d  . 

j  =  4  :  compatibility  condition  -^(s**  4-  =  Sx*«2  -  S*“3)  =  0.  (310) 


By  (3.9)  the  compatibility  condition  (3.10)  is  satisfied  identically.  The  compatibility  condition 
at  s  =  6  is  also  satisfied  identically  and  the  Kd\'  equation  is  thus  shown  to  be  Painleve. 

3.3  Macsyma  calculations 

We  want  to  carry  out  the  Painleve  test  for  the  Korteweg-de  Vries  equation.  Here  is  our  Macsyma 
session.  We  have  used  S  to  suppress  long  expressions  in  the  paper.  The  user  should  use  a  semicolon 
;  instead  to  see  the  results. 

(cl)  writefile(painlevetest)f 

Let  us  first  declare  the  dependencies  of  some  functions.  We  will  use 

(c2)  depends([f,g,pdlt,pdlx,pd3x,kdv,rec0,recl,rec2,rec3,rec4],[t,x])l 
We  define  f  to  be  the  following  finite  sum 

(c3)  f:g**-2*’sum(u[j](t,x)*g**j,  j,  0,  n); 

(d3) 

»  j=0 

Let  us  calculate  the  following  partial  derivatives  of  this  expression 
(c4)  pdlt:difr(f,t,l)8 
(c5)  pdlx:difr(f,x.l )S 
(c6)  pd3x:diff(f,x,3)S 

Now  we  add  these  terms  in  accordance  with  the  left-hand  side  of  the  KdV  equation 
(c7)  kdv  :  pdlt  -f  f*pdlx  -)-  b*pd3xt 

The  nonlinear  term  involves  the  product  f*(df/dx).  We  want  all  products  of  sums  to  be  converted 
in  nested  double  sums,  therefore  we  set 
(c8)  sumexpand:true8 
(c9)  gfactor(kdv)l 

Note  that  the  computer  returns  nested  sums  as  requested. 

There  is  a  common  denominator,  s',  which  we  remove  bv 
(clO)  %  •  g**5  I 

Let  us  put  factors,  such  as  powers  of  g,  inside  the  summations  by  setting 
(cll)  intosum(%)l 

We  have  to  choose  a  fixed  value  for  n,  for  our  purpose  n  =  6  suffices 


(cl2)  n:6  t 

(cl3)  ev(dll,suin)t 

We  calculate  the  coefficient  of  the  term  in  g°,  and  factor  the  result 
(cl4)  factor(ratcoef(dl3,g,0)); 

(dl4) 


(cl5)  cf:inpart(%,[^,2,3]); 
(dl5) 


(cl6)  rec0:dl4/cf; 
(dl6) 


-2gM<.«) 


uo(l,x)  +  l2fc(^)^ 

We  solve  for  Uo 

(cl7)  (.l)*part(%,2)f 
(cl8)  u[0](t,x)  :=  dl7l 
(cl9)  u[0](t,x)  ; 

(dl9) 

-I2Kg)’ 

Calculating  the  coefficent  of  g^,  substituting  for  uo  and  factoring  gives 
(c20)  factor(e v(ratcoef(dl3,g,l),diff)); 

(d20) 


306(^)»(uUf.x)-12fc^) 


(c21)  cf :  inpart(%,[l,2,3])S 
(c22)  reel  :  d20^f* 

(c23)  (-l)*part(%.2)l 
(c24)  ull](t,x)  :=  d23S 
(c25)  u[lj(t,x); 

(d28) 


Calculating  the  coefficient  of  g^.  substituting  Uq  and  uj  and  factoring  yields 
(c26)  factor(ev(ratcoef(dl3,g.2),difr))t 
(c27)  cf  :  inpart(%,[l,2,3])l 
(c28)  rec2  :  d26/cfS 

From  this  we  construct  U]  by  taking  parts,  etc. 

(c29)  inpart(%.allbut(4))t 
(c30)  (.l)*(%)/difr(g,x,l)*"2  S 


E'l:; 


(c31)  u[2](t,x)  :=  d30S 
(c32)  u[2](t,x); 

(d32) 
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Calculating  the  coefRcient  of  g^,  inserting  Uq,  and  V;  and  factoring  yields: 

(c33)  factor(ev(ratcoefr(dl3,g,3),difr)); 

(d33) 

Hi^^)  3(  ,  )  ^^dxdx^dz^  ^  dtdxdx^  didx^dx^^^dx 

This  expression  is  presumably  equivalent  to  the  recursion  relation  (3.9)  for  j  =  3  in  the  introduction. 

L''t  us  verify  this  by  doing  the  following  steps: 

(c34)  d33/(12*b*difr(g,x,l))# 

(c35)  expand((-l)*(%))8 
(c36)  combine(%)8 
(c37)  u[2](t,x)8 
(c38]  factor(part(d36,[3,4]))S 
(c39)  ratsubst(u2,d37,d38); 

(d39) 


(c40)  rec3  :  part(d36,[l,2,S])  +  d39; 

(d40) 

(c41)  u2  :  u[2](t,x)S 
(c42)  eqrecS  :  factor(ev(d40))S 
This  is  indeed  the  same  as  expressing  (d33)  after  division  by  (-12b) 

Let  us  construct  U3 

(c43)  (.l)*(%)»difr(g,x,l)**2l 
(c44)  iRpart(%,allbut(6))$ 

Note:  for  the  computer  the  first  term  in  (d43)  is  labelled  as  the  sixth  one! 

(c45)  u[3](t,x)  :=  (-l)*(d44)/difr(g,x,l)**4l 
(c46)  ui3](t,x); 

(d46) 

(S)* 

Finally,  we  verify  that  the  coefficient  of  g*  is  now  identically  zero,  after  all  the  substitutions  and 
simplifications  (compatibility  condition) 

(c47)  eqrec4  :  factor(ev(ratcoef(dl3,g,4),diff}); 

(d47)  0 


p  } Lin 


